#include <Functions.h>
#include <Tables.h>

//DFC for signal with fixed size, which can't be more then FURIE_MEMBERS_MAX
byte Furie(float* spectrum, bytep signal,  uint size)
{ 
	uint n = 0, k = 0;

	float _Complex[2] = {0, 0}; //Complex number, Complex[0] - real part, Complex[1] - imaginery part
	float _Complex_buf[2] = {0, 0}; 

	if(size > FURIE_MEMBERS_MAX)
		return 1;

	for(k = 0; k < size; k++)
	{
		for(n = 0; n < size; n++)
		{
			if(signal[n] == 0)
				continue;

			ComExp(n*k, size, _Complex_buf);		//Exponent (-2pi*n*k/N)
			ComMul(_Complex_buf, signal[n]);		//Multiply x[n] and exponent
			
			_Complex[0] += _Complex_buf[0];
			_Complex[1] += _Complex_buf[1];
		}
		spectrum[k] = ComModulus(_Complex);
	}
	
	return 0;
}

//Form windows on mel-scale
void MelPoints(float* window_Points)
{
	uint i = 0;
	float step = 0;

	window_Points[0] = 0;
	window_Points[WINDOW_COUNT] = GetMels(MAX_FREQUENCY);

	step = window_Points[WINDOW_COUNT]/(float)WINDOW_COUNT;

	for(i = 1; i < WINDOW_COUNT; i++)
		window_Points[i] = window_Points[i - 1] + step;
}

void Convolution(float* convolution, float* spectrum, float* windows_Points)
{
	uint i = 0, j = 0;
	
	float pre_point = 0;
	float post_point = 0;
	float central_point = 0;
	float delta_post = 0;
	float delta_pre = 0;

	float freq_j = 0;
	float freq_j_step = GetMels(FURIE_MEMBERS_MAX/(float)SAMPLE_RATE);

	float window = 0;

	for(i = 1; i <= WINDOW_COUNT; i++)
	{
		pre_point = windows_Points[i - 1];
		post_point = windows_Points[i + 1];
		central_point = windows_Points[i];
		delta_post = post_point - central_point;
		delta_pre = central_point - pre_point;

		freq_j = 0;
		convolution[i-1] = 0;

		for(j = 0; j < FURIE_MEMBERS_MAX; j++)
		{
			if(freq_j <= pre_point)
			{
				window = 0;

				freq_j += freq_j_step;
				continue;
			}

			if(freq_j > pre_point && freq_j <= central_point)
				window = (freq_j - pre_point)/delta_pre;

			if(freq_j > central_point && freq_j < post_point)
				window = (post_point - freq_j)/delta_post;

			if(freq_j >= post_point)
				break;

			convolution[i-1] += spectrum[j]*GetHerz(window);

			freq_j += freq_j_step;
		}
	}
}

void DCT(float* convolution, float* dct)
{
	uint i = 0, j = 0;
	
	for(i = 0; i < WINDOW_COUNT - 1; i++)
		convolution[i] = logf(convolution[i]);
	
	for(i = 0; i < DCT_MEMBERS; i++)
	{
		dct[i] = 0;

		for(j = 0; j < WINDOW_COUNT - 1; j++)
			dct[i] += convolution[j]*cosf((i+1)*(j+1-0.5)*pi/WINDOW_COUNT);
	}
}

float GetMels(float freq)
{
	float result = 1127.01048; //magical number with floating point, as we like
	result *= logf(1 + freq/700);
	
	return result;
}

float GetHerz(float freq)
{
	float result = 1127.01048; //magical number with floating point, as we like;
	result = 700*(expf(freq/result) - 1);

	return result;
}

//Exponent (-2pi*nom/denom)
void ComExp(uint num, uint denom, float* complex)
{
	complex[0] = Cos(num, denom)/65536.0;
	complex[1] = -Sin(num, denom)/65536.0;
}

//Multiply complex and real number
void ComMul(float* complex, int X)
{
	complex[0] = complex[0]*X;
	complex[1] = complex[1]*X;
}

//Modulus of complex number
float ComModulus(float* complex)
{
	return sqrtf(complex[0]*complex[0] + complex[1]*complex[1]);
}

//Cos(2pi*n/m)
short Cos(uint n, uint m)
{
	float nom = ceil(n*CScale/(float)m);

	if(nom <= CQuarter)
		return TCos[(int)nom];

	if(nom <= CHalfScale)
		return -(short)TCos[abs((int)(nom - CHalfScale))];

	if(nom <= CThreeQuarter)
		return -(short)TCos[(int)(nom - CHalfScale)];

	if(nom <= CScale)
		return TCos[abs((int)(nom - CScale))];
		
	return 0;
}

//Sin(2pi*n/m)
short Sin(uint n, uint m)
{
	int nom = (int)ceil(n*SScale/(float)m);

	if(nom <= SQuarter)
		return TSin[nom];

	if(nom <= SHalfScale)
		return TSin[abs(nom - SHalfScale)];

	if(nom <= SThreeQuarter)
		return -(short)TSin[nom - SHalfScale];

	if(nom <= SScale)
		return TSin[abs(nom - SScale)];

	return 0;
}

